\section{Simulation stuff}
In this section we summarize the set of equations that are used in
implementing the program to simulate the problem. The goal of the
program is to compute various physical quantities such as the energy,
heat capacity, magnetization, susceptibility that are associated with
our system. It is important to note that all these quantities are
finally reported on a ``per site'' basis.

The program considers $N_H$ spins in a 2-D triangular lattice with
$N_V$ layers that represent the time slices in the quantum
problem. One can effectively think of the system as 3-D lattice which
consists of many layers of the spins in a triangular lattice.

The total number of spins in our system is given by $N$ which is 
\begin{equation}
  \label{eq:11}
  N = N_V * N_H.
\end{equation}

The interaction strength in the vertical direction is given by
\begin{equation}
  \label{lambda}
  \lambda = - \frac 1 2 \ln \tanh(B \Delta \tau).
\end{equation}

The Trotter error goes as $ \Delta \tau^2 $ where 
\begin{equation}
  \label{dt}
  \Delta \tau = \frac {\beta} L = \frac 1 {T L}.
\end{equation}

The energy of the system can be thought of as consisting of two
separate components: (i) Energy of the spin-spin interaction denoted
$E_J$ and (ii) Energy of the interaction of the spins with the
magnetic field, $E_B$.

The initial values of $E_J$ and $E_B$ are given by
\begin{equation}
  \label{eq:21}
  E_J = - \frac 1 2 J \sum_{i \in N} S_i \eta_i^H
\end{equation}
and
\begin{equation}
  \label{eq:22}
  E_B = - B\sum_{i \in N_H, l \in N_V} e^{- 2 \lambda S_{il} S_{i (l+1)}} .
\end{equation}

The action of our system i.e. the 3-D spin lattice is
\begin{equation}
  \label{eq:23}
  S = - J \Delta \tau \sum_i S_i \eta^H_i - \lambda \sum_i S_i \eta^V_i 
=  \sum_i S_i ( - J \Delta \tau \eta^H_i - \lambda  \eta^V_i),
\end{equation}
where the variables $\eta^H_i$ and $\eta^V_i$ denote the sum of the
horizontal and vertical neighbors of the spin $i$ respectively.For the
above action, (according to the Metropolis algorithm, ) we obtain the
probability of flipping a spin $i$ is
\begin{equation}
  \label{eq:24}
  P_i = \frac {e^{-\Delta S}} {1 +e^{-\Delta S} }
\end{equation}
Where 
\begin{equation}
  \label{eq:dS}
  \Delta S = 2 S_i ( J \Delta \tau \eta^H_i + \lambda  \eta^V_i)
\end{equation}

If the spin flip is flipped the energies associated with both the
spin-spin interaction ($E^{\phantom \dagger}_J$) and the spin-field interaction ($E^{\phantom \dagger}_B$) will
change. Expressions for the change in these quantities in terms of the
value of the spin before the flip takes place is,
\begin{equation}
  \label{eq:25}
  \Delta E_J(i) = 2 J S_i \eta^H_i,
\end{equation}

\begin{gather}
  \Delta E_B (i) = - 2 B (\sinh( 2 \lambda S_{il} S_{i (l+1)}) + \sinh(2 \lambda
  S_{il} S_{i (l-1)}))=\\ = -2 B \eta^V_i \sinh(2 \lambda S_{il} S_{i(l+1)}).
\end{gather}

\subsection{Tests to check program}

1. If $B = 0$, $\lambda = 0$, $N_V = 1$ It should represent classical 1 layer problem\\
2. Initally we calculate $E_J$. Then each step we add some $\Delta E_J$. After $K$ steps exect expression for $E_J$ and that one we got from $E_J + \Delta E_{J1}\ldots \Delta E_{JK}$ should be equal. The same for $E_B$. $|\Delta E_B| = 0$ or const on each step\\
3. If $J = 0$ We can calculate exect value of average energy.

\begin{gather}
  Z = e^{-\beta B} + e^{\beta B}\\
E = \frac B Z (e^{-\beta B} - e^{\beta B})=\\
= B \frac {e^{-\beta B} - e^{\beta B}} {e^{-\beta B} + e^{\beta B}}\\
= - B \tanh (\beta B)
\end{gather}